A Quick Start Tutorial

This tutorial offers a quick-start guide for using BaySpec to fit a spectral model to gamma-ray data. It can be broadly divided into the following three sections:

• Data: The spectra from Fermi/GBM’s NaI detector and BGO detector.

• Model: A simple cutoff power-law function.

• Fitting: Bayesian inference implemented using multinest.

import numpy as np
from bayspec.model.local import *
from bayspec import DataUnit, Data, Infer, Plot

1. Load spectra and response data.

nai = DataUnit(
    src='./ME/me.src', 
    bkg='./ME/me.bkg', 
    rsp='./ME/me.rsp', 
    notc=[8, 900], 
    stat='pgstat', 
    grpg={'min_sigma': 3, 'max_bin': 10})

bgo = DataUnit(
    src='./HE/he.src', 
    bkg='./HE/he.bkg', 
    rsp='./HE/he.rsp', 
    notc=[300, 38000], 
    stat='pgstat', 
    grpg={'min_sigma': 3, 'max_bin': 10})

data = Data([('nai', nai), 
             ('bgo', bgo)])

2.Define spectral model.

model = cpl()

3.Run bayesian infer.

infer = Infer([(data, model)])
print(infer)

╒════════╤══════════════╤═════════════╤═════════════╤═════════╕
│  cfg#  │  Expression  │  Component  │  Parameter  │  Value  │
╞════════╪══════════════╪═════════════╪═════════════╪═════════╡
│   1    │     cpl      │     cpl     │  redshift   │    0    │
╘════════╧══════════════╧═════════════╧═════════════╧═════════╛
╒════════╤══════════════╤═════════════╤═════════════╤═════════╤═════════════╕
│  par#  │  Expression  │  Component  │  Parameter  │  Value  │    Prior    │
╞════════╪══════════════╪═════════════╪═════════════╪═════════╪═════════════╡
│   1*   │     cpl      │     cpl     │  $\alpha$   │   -1    │ unif(-8, 4) │
├────────┼──────────────┼─────────────┼─────────────┼─────────┼─────────────┤
│   2*   │     cpl      │     cpl     │ log$E_{c}$  │    2    │ unif(0, 4)  │
├────────┼──────────────┼─────────────┼─────────────┼─────────┼─────────────┤
│   3*   │     cpl      │     cpl     │   log$A$    │   -1    │ unif(-6, 5) │
╘════════╧══════════════╧═════════════╧═════════════╧═════════╧═════════════╛

post = infer.multinest(nlive=300, resume=True, savepath='./multinest')

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  300
 dimensionality =    3
 resuming from previous job
 *****************************************************
 Starting MultiNest
Acceptance Rate:                        0.689121
Replacements:                               6271
Total Samples:                              9100
Nested Sampling ln(Z):               -229.513355
Importance Nested Sampling ln(Z):    -229.263794 +/-  0.024273
 ln(ev)=  -229.21088552733494      +/-  0.24113203804165545     
 Total Likelihood Evaluations:         9100
 Sampling finished. Exiting MultiNest
  analysing data from ./multinest/1-.txt

print(post)

╒════════╤══════════════╤═════════════╤═════════════╤════════╤══════════╤════════╤══════════════════╕
│  par#  │  Expression  │  Component  │  Parameter  │  Mean  │  Median  │  Best  │    1sigma CI     │
╞════════╪══════════════╪═════════════╪═════════════╪════════╪══════════╪════════╪══════════════════╡
│   1    │     cpl      │     cpl     │  $\alpha$   │ -1.563 │  -1.563  │ -1.561 │ [-1.572, -1.552] │
├────────┼──────────────┼─────────────┼─────────────┼────────┼──────────┼────────┼──────────────────┤
│   2    │     cpl      │     cpl     │ log$E_{c}$  │  2.69  │   2.69   │ 2.688  │  [2.673, 2.707]  │
├────────┼──────────────┼─────────────┼─────────────┼────────┼──────────┼────────┼──────────────────┤
│   3    │     cpl      │     cpl     │   log$A$    │ -0.771 │  -0.771  │ -0.77  │ [-0.778, -0.765] │
╘════════╧══════════════╧═════════════╧═════════════╧════════╧══════════╧════════╧══════════════════╛
╒════════╤═════════╤═════════════╤════════════╤════════╕
│  Data  │  Model  │  Statistic  │   Value    │  Bins  │
╞════════╪═════════╪═════════════╪════════════╪════════╡
│  nai   │   cpl   │   pgstat    │   388.00   │  106   │
├────────┼─────────┼─────────────┼────────────┼────────┤
│  bgo   │   cpl   │   pgstat    │   32.62    │   26   │
├────────┼─────────┼─────────────┼────────────┼────────┤
│ Total  │  Total  │  stat/dof   │ 420.62/129 │  132   │
╘════════╧═════════╧═════════════╧════════════╧════════╛
╒════════╤════════╤════════╤═════════╕
│  AIC   │  AICc  │  BIC   │   lnZ   │
╞════════╪════════╪════════╪═════════╡
│ 426.62 │ 426.81 │ 435.27 │ -229.26 │
╘════════╧════════╧════════╧═════════╛

fig = Plot.infer_ctsspec(post, style='CE', show=True)

fig = Plot.post_corner(post, show=True)

earr = np.logspace(np.log10(0.5), 3, 100)

modelplot = Plot.model(ploter='plotly', style='vFv', CI=True)
fig = modelplot.add_model(model, E=earr, show=True)